/** @file essl_impl_impl.h
 * @brief Implementation of essl_impl-enabled FFT engines
 *
 */
#ifndef QMCPLUSPLUS_FFT_ESSL_IMPL_H
#define QMCPLUSPLUS_FFT_ESSL_IMPL_H
#include <essl.h>
#include <iostream>
#include <vector>

template<class T>
struct fft_essl_plan_type
{
  T* aux1;
  T* aux2;

  fft_essl_plan_type():aux1(0),aux2(0){}

  ~fft_essl_plan_type()
  {
    clean();
  }

  void clean()
  {
    if(aux1) delete [] aux1; aux1=0;
    if(aux2) delete [] aux2; aux2=0;
  }

  void resize(int n1, int n2)
  {
    aux1=new double[n1];
    aux2=new double[n2];
  }
};

/** base class for FFT engines using essl library
 */
struct fft_essl_base
{
  static const int FFTESSL_FORWARD=1;
  static const int FFTESSL_BACKWARD=-1;

  typedef fft_essl_plan_type<double> fft_plan_type;

  int naux1;
  int naux2;
  fft_plan_type forward_plan;
  fft_plan_type backward_plan;

  ///default constructor
  fft_essl_base() 
  { }

  ///virtual destructor to clean up the plans
  virtual ~fft_essl_base()
  {
    clean();
  }

  void clean()
  {
    forward_plan.clean();
    backward_plan.clean();
  }
};


/** specialization for single-fft only with essl-FFT on complex
 */
template<>
struct FFT_DFT<double,std::complex<double>,false,FFTESSL_ENG> 
: public fft_base, public fft_essl_base
{
  typedef std::complex<double> value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,1,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,1);
      fft_essl_base::clean();

      naux1=(fft_size>2048)?static_cast<int>(20000.+2.28*dims)+1:20000;
      naux2=naux1;

      forward_plan.resize(naux1,naux2);
      backward_plan.resize(naux1,naux2);

      dcft(1,in,1,1,in,1,1,dims,1,FFTESSL_FORWARD,1.0
          ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);

      dcft(1,in,1,1,in,1,1,dims,1,FFTESSL_BACKWARD,1.0
          ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
    }
  }

  inline void transform_forward_atom(value_type* in)
  {
    dcft(0,in,1,1,in,1,1,fft_size,1,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void transform_backward_atom(value_type* in)
  {
    dcft(0,in,1,1,in,1,1,fft_size,1,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1 ,backward_plan.aux2,naux2);
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) 
      transform_forward_atom(offset);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_size) 
      transform_backward_atom(offset);
  }
};

/** specialization for single-fft only with essl-FFT on double
 */
template<>
struct FFT_DFT<double,double,false,FFTESSL_ENG> 
: public fft_base, public fft_essl_base
{
  typedef double value_type;
  int fft_mem_offset;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,1,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_mem_offset=dims+2;
      fft_base::resize(dims,1);
      fft_essl_base::clean();

      naux1=(fft_size>4096)?static_cast<int>(20000.+1.64*dims)+1:25000;
      naux2=naux1;
      forward_plan.resize(naux1,naux2);
      backward_plan.resize(naux1,naux2);

      drcft(1,in,1,in,1,fft_size,1,FFTESSL_FORWARD,1.0
          ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);

      dcrft(1,in,1,in,1,fft_size,1,FFTESSL_BACKWARD,1.0
          ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
    }
  }

  inline void transform_forward_atom(value_type* in)
  {
    drcft(0,in,1,in,1,fft_size,1,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void transform_backward_atom(value_type* in)
  {
    dcrft(0,in,1,in,1,fft_size,1,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_mem_offset) 
      transform_forward_atom(offset);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    value_type* offset=in;
    for(int i=0; i<howmany; ++i,offset+=fft_mem_offset) 
      transform_backward_atom(offset);
  }
};

/** specialization for multiple-fft with essl-FFT on complex
 */
template<>
struct FFT_DFT<double,std::complex<double>,true,FFTESSL_ENG> 
: public fft_base, public fft_essl_base
{
  typedef std::complex<double> value_type;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,howmany,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      fft_base::resize(dims,howmany);
      fft_essl_base::clean();

      naux1=(fft_size>2048)?static_cast<int>(20000.+2.28*dims)+1:20000;
      naux2=std::max(naux1,(2*dims+256)*std::min(64,howmany));

      forward_plan.resize(naux1,naux2);
      backward_plan.resize(naux1,naux2);

      dcft(1,in,1,dims,in,1,dims,dims,howmany,FFTESSL_FORWARD,1.0
          ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);

      dcft(1,in,1,dims,in,1,dims,dims,howmany,FFTESSL_BACKWARD,1.0
          ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    dcft(0,in,1,fft_size,in,1,fft_size,fft_size,howmany,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    dcft(0,in,1,fft_size,in,1,fft_size,fft_size,howmany,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }
};

/** specialization for single-fft only with essl-FFT on double
 */
template<>
struct FFT_DFT<double,double,true,FFTESSL_ENG> 
: public fft_base, public fft_essl_base
{
  typedef double value_type;
  int forward_offset;
  int backward_offset;

  FFT_DFT(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    init(dims,in,howmany,uflag);
  }

  void init(int dims, value_type* in, int howmany=1, unsigned uflag=FFTW_MEASURE)
  {
    if(dims != fft_size)
    {
      forward_offset=dims+2;
      backward_offset=dims/2+1;
      fft_base::resize(dims,howmany);
      fft_essl_base::clean();

      naux1=(fft_size>4096)?static_cast<int>(20000.+1.64*dims)+1:25000;
      naux2=naux1;
      forward_plan.resize(naux1,naux2);
      backward_plan.resize(naux1,naux2);

      drcft(1,in,forward_offset,in,backward_offset,fft_size,howmany,FFTESSL_FORWARD,1.0
          ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);

      dcrft(1,in,backward_offset,in,forward_offset,fft_size,howmany,FFTESSL_BACKWARD,1.0
          ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
    }
  }

  inline void transform_forward(value_type* in, int howmany=1)
  {
    drcft(0,in,forward_offset,in,backward_offset,fft_size,howmany,FFTESSL_FORWARD,1.0
        ,forward_plan.aux1,naux1 ,forward_plan.aux2,naux2);
  }

  inline void transform_backward(value_type* in, int howmany=1)
  {
    dcrft(0,in,backward_offset,in,forward_offset,fft_size,howmany,FFTESSL_BACKWARD,1.0
        ,backward_plan.aux1,naux1,backward_plan.aux2,naux2);
  }
};
#endif
